#
# PlotIRTResult.R
#
# Replication file for
# Hill, Seth J. and Gregory A. Huber. 2018. ``On The Meaning of Survey Reports of Roll Call `Votes','' American Journal of Political Science.
#

rm(list=ls())
if (.Platform$OS == "windows") {
  # Set working directory in location of script.
  .doit <- function() { # only works with R.exe; trap errors if using Rscript.exe
  frame_files <- lapply(sys.frames(), function(x) x$ofile);frame_files <- Filter(Negate(is.null), frame_files)
  PATH <- dirname(frame_files[[length(frame_files)]]); setwd(PATH) ; rm(PATH,frame_files)
  }
  try(.doit(),silent=T)
}

library(foreign)
library(pscl)
library(data.table)
library(Hmisc) # for wtd.quantile.
library(RColorBrewer)
palette(brewer.pal(n=9,name="Set1")[c(1:5,7:9)]) 
par(cex.main=1.2)

iters <- -1 # set to -1 to use all iterations
if (iters != -1) warning("Not using all iterations to summarize.")
.anArrow <- function(x0,y0,x1=x0,y1,length=.2,label=NULL,...) {
  # Make an arrow on a plot with optional text.
  arrows(x0,y0,x1,y1,length=length,...)
  if (!is.null(label)) {
    text(x0,y0,label=label,...)
  }
}

#
# Load results from IrtScaleVotes.R.
#
f <- ("RawData/IRTModelResults-Screened.RData")
load(file=f)

#
# Turn results into data.tables.
#
# Nb: ideal() result object has list element
# `x' that, per documentation, is "a three-dimensional 
# array containing the MCMC output with respect to the
# the ideal point of each legislator in each dimension.
# The three-dimensional array is in 
# iteration-legislator-dimension order. The iterations
# run from burnin to maxiter, at an interval of thin."
#
cat("Creating long data.tables...\n");flush.console()
.makeLong <- function(x,iters=-1) {
  # Create a long data.table with first dimension
  # ideal point estimates from ideal result x.
  # Arguments.
  #  x -- the x object from and ideal() result.
  #  iters -- number of iterations to store; -1 stores all.
  res <- NULL
  if (iters == -1) { iters <- dim(x)[1] }
  for (i in 1:dim(x)[2]) { # loop over legislators
    res <- rbindlist(list(res,data.table(iter=seq_len(iters),id=rep(dimnames(x)[[2]][i],iters),x=x[1:iters,i,1])))
  }
  res
}
long.sen <- .makeLong(sen.id$id$x,iters)
long.std <- .makeLong(std.id$id$x,iters)
long.inf <- .makeLong(inf.id$id$x,iters)
long.sen.std <- .makeLong(sen.std.id$id$x,iters)
long.sen.inf <- .makeLong(sen.inf.id$id$x,iters)
# Join variables weight, issen, and sen.pid to 
# long data.tables.
DT <- data.table(dat) # issen and sen.pid
DT[,id := as.character(V1)]
# Merge weight estimated in RakeToPew.R.
wts <- data.table(read.dta("RawData/RakedPewWeights.dta"))
DT <- merge(DT,wts,by="V1",all.x=T,all.y=F)
DT[is.na(weight), weight := 1] # make Sen weights 1
# Merge to long tables.
long.sen <- merge(long.sen,DT[,c("id","is.sen","sen.pid","weight"),with=F],by="id",all.x=T,all.y=F)
long.std <- merge(long.std,DT[,c("id","weight"),with=F],by="id",all.x=T,all.y=F)
long.inf <- merge(long.inf,DT[,c("id","weight"),with=F],by="id",all.x=T,all.y=F)
long.sen.std <- merge(long.sen.std,DT[,c("id","is.sen","sen.pid","weight"),with=F],by="id",all.x=T,all.y=F)
long.sen.inf <- merge(long.sen.inf,DT[,c("id","is.sen","sen.pid","weight"),with=F],by="id",all.x=T,all.y=F)

#
# Rotation to conservative = positive.
#
cat("Checking on rotation...\n");flush.console()

# Determine if scaling put conservative in positive
# direction by checking discrimination param for
# repeal ACA.
.checkACA <- function(id) {
  # Return sign of repeal ACA discrim param.
  sign(id$betabar["bill.rpl.aca",1])
}
.rotateConservative <- function(long,id) {
  # Check if ideal points need to be rotated to
  # make conservative positive.
  sg <- .checkACA(id) # should be positive
  if (sg == -1) {
    cat("Flipping ideal points so positive=conservative.\n")
    long[,x := -1*x]
  }
}
.rotateConservative(long.sen,sen.id$id)
.rotateConservative(long.std,std.id$id)
.rotateConservative(long.inf,inf.id$id)
.rotateConservative(long.sen.std,sen.std.id$id)
.rotateConservative(long.sen.inf,sen.inf.id$id)

#
# Plot of posterior median densities for
# Senators and respondents, separately by
# condition.
#
cat("Making posterior densities plots...\n");flush.console()
.pDensPlot <- function(long,med=T,add=F,lcol="black",...) {
  # Plot the density of posterior samples.
  # Arguments.
  #  long -- data.table of posterior samples.
  #  med -- plot density of posterior medians (when T), 
  #         or of all posterior draws (when F).
  #  add -- add lines to current plot window when T.
  #  lcol -- line color.
  #  ... -- other arguments passed to plot().
  if (med) {
    meds <- long[,median(x,na.rm=T),by=c("id","weight")]
  } else {
    meds <- copy(long[,c("id","weight","x"),with=F])
    setnames(meds,gsub("x","V1",names(meds)))
  }
  meds[,weight := weight/sum(weight)]
  if (add) { # Add lines
    dens <- meds[,density(x=V1,weights=weight,bw="nrd0",adjust=1.2)]
    # Rescale max of dens to par('usr')[4]
    dens$y <- par('usr')[4]*dens$y/max(dens$y)
    lines(dens,lwd=3,col=lcol,...)
  } else { # New plot
    plot(meds[,density(x=V1,weights=weight,bw="nrd0",adjust=1.2)],type='l',lwd=3,col=lcol,ann=F,axes=F,...)
    axis(1);title(xlab="Ideal point",line=2)
  }
}

#
# Posterior Senator and respondent quantiles, 
# separately by respondent condition.
#
cat("Plotting posterior quantiles...\n");flush.console()

.calcQuants <- function(long,quant=.5,byv="is.sen") {
  # Return a data.table with the quantile of the ideal
  #  point distribution on each posterior iteration.
  # Arguments.
  #  long -- data.table with posterior draws
  #  quant -- desired quantile
  #  byv -- by variable.

  library(Hmisc) # for wtd.quantile.
  # Quantile of distribution on each iter by byv.
  qn <- long[,list(qnt=wtd.quantile(x,probs=quant,weights=weight)),by=c(byv,"iter")]
  qn
}
#d <- .calcQuants(long.sen.std)

# Calculate multiple quantiles, both for
# standard and info scalings of respondents vs Senators.
quants <- c(.1,.25,.41,.5,.6,.75,.9)
res <- NULL
cat("Calculating quantiles...\n");flush.console()
for (qn in quants) {
  cat(" Quantile",qn,"\n");flush.console()
  # Quantile for standard condition.
  d <- .calcQuants(long.sen.std,quant=qn)
  d[,cond := "Standard"];d[,pr := qn]
  res <- rbindlist(list(res,d))
  # Distance for info condition.
  d <- .calcQuants(long.sen.inf,quant=qn)
  d[,cond := "Info"];d[,pr := qn]
  res <- rbindlist(list(res,d))
}

# Calculate posterior median distance with 95
# pct credible interval.
qns <- res[,list(med=median(qnt),lb=quantile(qnt,.025),ub=quantile(qnt,.975)),by=c("pr","cond","is.sen")]

# Plot.
pdf("figures/FigureA02IRTQuantiles-Screened.pdf",width=11,height=6)
par.old <- par(cex.axis=1.2,cex.lab=1.2,mar=c(4.1,4.1,2.1,1.1))
# Standard condition.
.quantPlot <- function(qns,quants) {
  # Plot ideal point and 95% CI for quantiles.
  qns[,plot(x=pr+ifelse(is.sen==1,-0.01,0.01),y=med,pch=ifelse(is.sen==1,19,15),cex=2,ylim=range(c(lb,ub)),ann=F,axes=F,xlim=c(0.05,1))]
  legend('topleft',legend=c("Senators","Respondents"),pch=c(19,15),cex=1.2)
  axis(1,at=quants);axis(2,las=2)
  abline(h=0,lty=2,lwd=2)
  qns[,segments(x0=pr+ifelse(is.sen==1,-0.01,0.01),y0=lb,y1=ub,lwd=2)]
  title(ylab="Posterior quantile")
  title(xlab="Quantile of ideal point distribution",line=3)
  # Up arrow.
  .anArrow(x0=.99,y0=0.05,y1=1.4,lwd=2)
  text(x=.98,y=1.1,pos=2,srt=90,labels="More conservative",cex=1.2)
  # Down arrow.
  .anArrow(x0=.99,y0=-0.05,y1=-1.4,lwd=2)
  text(x=.98,y=-.15,pos=2,srt=90,labels="More liberal",cex=1.2)
}
.quantPlot(qns[cond=="Standard",],quants)
title(main="Joint scaling, standard condition",line=1)
# Info condition.
.quantPlot(qns[cond=="Info",],quants)
title(main="Joint scaling, informational condition",line=1)
par(par.old)
dev.off()

# Write quantiles to csv for reference in writing.
write.csv(qns,"PosteriorQuantilesByCondAndSen-Screened",row.names=F)




# Unscreened


#
# Load results from IrtScaleVotes.R.
#
f <- ("RawData/IRTModelResults-Unscreened.RData")
load(file=f)

#
# Turn results into data.tables.
#
# Nb: ideal() result object has list element
# `x' that, per documentation, is "a three-dimensional 
# array containing the MCMC output with respect to the
# the ideal point of each legislator in each dimension.
# The three-dimensional array is in 
# iteration-legislator-dimension order. The iterations
# run from burnin to maxiter, at an interval of thin."
#
cat("Creating long data.tables...\n");flush.console()
.makeLong <- function(x,iters=-1) {
  # Create a long data.table with first dimension
  # ideal point estimates from ideal result x.
  # Arguments.
  #  x -- the x object from and ideal() result.
  #  iters -- number of iterations to store; -1 stores all.
  res <- NULL
  if (iters == -1) { iters <- dim(x)[1] }
  for (i in 1:dim(x)[2]) { # loop over legislators
    res <- rbindlist(list(res,data.table(iter=seq_len(iters),id=rep(dimnames(x)[[2]][i],iters),x=x[1:iters,i,1])))
  }
  res
}
long.sen <- .makeLong(sen.id$id$x,iters)
long.std <- .makeLong(std.id$id$x,iters)
long.inf <- .makeLong(inf.id$id$x,iters)
long.sen.std <- .makeLong(sen.std.id$id$x,iters)
long.sen.inf <- .makeLong(sen.inf.id$id$x,iters)
# Join variables weight, issen, and sen.pid to 
# long data.tables.
DT <- data.table(dat) # issen and sen.pid
DT[,id := as.character(V1)]
# Merge weight estimated in RakeToPew.R.
wts <- data.table(read.dta("RawData/RakedPewWeights.dta"))
DT <- merge(DT,wts,by="V1",all.x=T,all.y=F)
DT[is.na(weight), weight := 1] # make Sen weights 1
# Merge to long tables.
long.sen <- merge(long.sen,DT[,c("id","is.sen","sen.pid","weight"),with=F],by="id",all.x=T,all.y=F)
long.std <- merge(long.std,DT[,c("id","weight"),with=F],by="id",all.x=T,all.y=F)
long.inf <- merge(long.inf,DT[,c("id","weight"),with=F],by="id",all.x=T,all.y=F)
long.sen.std <- merge(long.sen.std,DT[,c("id","is.sen","sen.pid","weight"),with=F],by="id",all.x=T,all.y=F)
long.sen.inf <- merge(long.sen.inf,DT[,c("id","is.sen","sen.pid","weight"),with=F],by="id",all.x=T,all.y=F)

#
# Rotation to conservative = positive.
#
cat("Checking on rotation...\n");flush.console()

# Determine if scaling put conservative in positive
# direction by checking discrimination param for
# repeal ACA.
.checkACA <- function(id) {
  # Return sign of repeal ACA discrim param.
  sign(id$betabar["bill.rpl.aca",1])
}
.rotateConservative <- function(long,id) {
  # Check if ideal points need to be rotated to
  # make conservative positive.
  sg <- .checkACA(id) # should be positive
  if (sg == -1) {
    cat("Flipping ideal points so positive=conservative.\n")
    long[,x := -1*x]
  }
}
.rotateConservative(long.sen,sen.id$id)
.rotateConservative(long.std,std.id$id)
.rotateConservative(long.inf,inf.id$id)
.rotateConservative(long.sen.std,sen.std.id$id)
.rotateConservative(long.sen.inf,sen.inf.id$id)

#
# Plot of posterior median densities for
# Senators and respondents, separately by
# condition.
#
cat("Making posterior densities plots...\n");flush.console()
.pDensPlot <- function(long,med=T,add=F,lcol="black",...) {
  # Plot the density of posterior samples.
  # Arguments.
  #  long -- data.table of posterior samples.
  #  med -- plot density of posterior medians (when T), 
  #         or of all posterior draws (when F).
  #  add -- add lines to current plot window when T.
  #  lcol -- line color.
  #  ... -- other arguments passed to plot().
  if (med) {
    meds <- long[,median(x,na.rm=T),by=c("id","weight")]
  } else {
    meds <- copy(long[,c("id","weight","x"),with=F])
    setnames(meds,gsub("x","V1",names(meds)))
  }
  meds[,weight := weight/sum(weight)]
  if (add) { # Add lines
    dens <- meds[,density(x=V1,weights=weight,bw="nrd0",adjust=1.2)]
    # Rescale max of dens to par('usr')[4]
    dens$y <- par('usr')[4]*dens$y/max(dens$y)
    lines(dens,lwd=3,col=lcol,...)
  } else { # New plot
    plot(meds[,density(x=V1,weights=weight,bw="nrd0",adjust=1.2)],type='l',lwd=3,col=lcol,ann=F,axes=F,...)
    axis(1);title(xlab="Ideal point",line=2)
  }
}

#
# Posterior Senator and respondent quantiles, 
# separately by respondent condition.
#
cat("Plotting posterior quantiles...\n");flush.console()

.calcQuants <- function(long,quant=.5,byv="is.sen") {
  # Return a data.table with the quantile of the ideal
  #  point distribution on each posterior iteration.
  # Arguments.
  #  long -- data.table with posterior draws
  #  quant -- desired quantile
  #  byv -- by variable.

  library(Hmisc) # for wtd.quantile.
  # Quantile of distribution on each iter by byv.
  qn <- long[,list(qnt=wtd.quantile(x,probs=quant,weights=weight)),by=c(byv,"iter")]
  qn
}
#d <- .calcQuants(long.sen.std)

# Calculate multiple quantiles, both for
# standard and info scalings of respondents vs Senators.
quants <- c(.1,.25,.41,.5,.6,.75,.9)
res <- NULL
cat("Calculating quantiles...\n");flush.console()
for (qn in quants) {
  cat(" Quantile",qn,"\n");flush.console()
  # Quantile for standard condition.
  d <- .calcQuants(long.sen.std,quant=qn)
  d[,cond := "Standard"];d[,pr := qn]
  res <- rbindlist(list(res,d))
  # Distance for info condition.
  d <- .calcQuants(long.sen.inf,quant=qn)
  d[,cond := "Info"];d[,pr := qn]
  res <- rbindlist(list(res,d))
}

# Calculate posterior median distance with 95
# pct credible interval.
qns <- res[,list(med=median(qnt),lb=quantile(qnt,.025),ub=quantile(qnt,.975)),by=c("pr","cond","is.sen")]

# Plot.
pdf("figures/FigureA06IRTQuantiles-Unscreened.pdf",width=11,height=6)
par.old <- par(cex.axis=1.2,cex.lab=1.2,mar=c(4.1,4.1,2.1,1.1))
# Standard condition.
.quantPlot <- function(qns,quants) {
  # Plot ideal point and 95% CI for quantiles.
  qns[,plot(x=pr+ifelse(is.sen==1,-0.01,0.01),y=med,pch=ifelse(is.sen==1,19,15),cex=2,ylim=range(c(lb,ub)),ann=F,axes=F,xlim=c(0.05,1))]
  legend('topleft',legend=c("Senators","Respondents"),pch=c(19,15),cex=1.2)
  axis(1,at=quants);axis(2,las=2)
  abline(h=0,lty=2,lwd=2)
  qns[,segments(x0=pr+ifelse(is.sen==1,-0.01,0.01),y0=lb,y1=ub,lwd=2)]
  title(ylab="Posterior quantile")
  title(xlab="Quantile of ideal point distribution",line=3)
  # Up arrow.
  .anArrow(x0=.99,y0=0.05,y1=1.4,lwd=2)
  text(x=.98,y=1.1,pos=2,srt=90,labels="More conservative",cex=1.2)
  # Down arrow.
  .anArrow(x0=.99,y0=-0.05,y1=-1.4,lwd=2)
  text(x=.98,y=-.15,pos=2,srt=90,labels="More liberal",cex=1.2)
}
.quantPlot(qns[cond=="Standard",],quants)
title(main="Joint scaling, standard condition",line=1)
# Info condition.
.quantPlot(qns[cond=="Info",],quants)
title(main="Joint scaling, informational condition",line=1)
par(par.old)
dev.off()

# Write quantiles to csv for reference in writing.
write.csv(qns,"PosteriorQuantilesByCondAndSen-Unscreened",row.names=F)


